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Abstract. We present a new model of giant planet formation that extends the core-accretion model of Pollack et 
0.1- ( 1996J to include migration, disc evolution and gap formation. We show that taking into account these effects 
can lead to a much more rapid formation of giant planets, making it compatible with the typical disc lifetimes 
inferred from observations of young circumstellar discs. This speed up is due to the fact that migration prevents 
the severe depletion of the feeding zone as observed in in situ calculations. Hence, the growing planet is never 
isolated and it can reach cross-over mass on a much shorter timescale. To illustrate the range of planets that can 
form in our model, we describe a set of simulations in which we have varied some of the initial parameters and 
compare the final masses and semi-major axes with those inferred from observed extra-solar planets. 
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1. Introduction 

The standard giant planet formation scenario is the so- 
called core-accretion model. In this model, a solid core 
is formed first by accretion of solid planetesimals which 
themselves were formed by sedimentation and coagulation 
of small dust grains (Wetherill & Steward 1989, Lissauer 
1993). As the core grows, it eventually becomes massive 
enough to gravitationally bind some of the nebular gas 
thus, surrounding itself by a tenuous envelope. The subse- 
quent evolution of this core/envelope structure has been 
studied first by Perri & Cameron (1974) and subsequently 
in great detail by Bodenheimer & Pollack (1986 hereafter 
referred to as BP86) and Pollack et al. (1996 hereafter 
referred to as P96). These authors have found that both 
the solid core and the gaseous envelope subsequently grow 
in mass, the envelope remaining in quasi-static and ther- 
mal equilibrium. During this phase, the energy radiated 
by the gas is supplied by the energy released from the 
accretion of planetesimals. As the planet reaches a high 
enough mass (of the order of 20 — SOM^ at 5 AU, but de- 
pending on different physical parameters such as the solid 
accretion rate) , radiative losses become so large that they 
can no longer be offset by planetesimal accretion alone and 
the envelope starts to contract. At this stage, the mass is 
nearly equally distributed between the mass of accreted 
gas, and the mass of accreted planetesimals. The contrac- 
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tion of the envelope increases the gas accretion rate which 
in turn increases the energy losses through radiation and 
the process runs away very rapidly building up a massive 
envelope. 

Detailed numerical calculations in particular by P96 
have shown that this core-accretion formation model can 
be subdivided in three distinct phases. During phase 1, 
the core is formed by accretion of planetesimals located 
inside the growing planet's feeding zone which extends to 
a few Hills radii. The accretion rate of solids during this 
phase is typically of the order of lO~^M0/yr while the 
gas accretion rate is several orders of magnitude lower. 
Phase 1 ends when the feeding zone of the planet becomes 
severely depleted, which generally occurs before the planet 
has reached a mass high enough to accrete a large amount 
of gas. During phase 2, the mass increase is essentially 
due to the slow accretion of gas in the envelope. Note that 
by increasing the envelope mass, the feeding zone of the 
planet also increases which in turns allows the accretion of 
more planetesimals. Both accretion rates (solid and gas) 
turn out to be relatively constant during this phase with 
the gas accretion rate exceeding by a fraction of an or- 
der of magnitude the solid accretion rate. This phase lasts 
until the planet has reached the cross-over mass (mass 
of accreted planetesimals equal mass of accreted gas), at 
which time the system enters phase 3. At this point, evo- 
lution proceeds extremely fast by runaway gas accretion 
as the envelope is no longer able to maintain quasi-static 
equilibrium. The mass of the planet increases correspond- 
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ingly. The timescale for the formation of a giant planet in 
this core-accretion scenario is almost entirely determined 
by phase 2, which, as shown by P96, turns out to be an 
extremely sensitive function of the disc surface density. 
P96 found for their preferred model for the formation of 
Jupiter (model Jl, the surface density being equal to ^ 4 
times the value in the minimum mass solar nebula) a for- 
mation timescale close to 8 Myr while reducing the surface 
density to 75% of this value leads to a formation time of 
nearly 50 Myr. 

The core-accretion scenario has been motivated by the 
apparent existence of a solid core, long estimated to be of 
the order of the mass of accreted planetesimals when the 
planet reaches the cross-over mass, in all the giant planets 
of the solar system. This scenario has also been supported 
by the enrichment in heavy elements (heavier than He) of 
both Jupiter and Saturn compared to solar value, deduced 
from interior structure models and various remote/in situ 
measurements (radius, mass, surface abundance, gravita- 
tional moments, see Guillot et al. 2004). 

However, the major difficulty affecting the core- 
accretion scenario which has been pointed out repeat- 
edly is related to the timescale required to form a giant 
planet. Based on astronomical observations, protoplane- 
tary discs are believed to transport mass inward at a rate 
of IQ-^^^^Mg/yr (Hartmann et al. 1998) while their total 
mass has been estimated to lie between 10~^ and 1O~'^M0 
(Beckwith & Sargent 1996). From these numbers as well 
as from the observations that about half the stars in young 
clusters loose their discs within 3 Myr (Haisch et al. 2001) 
one infers a typical circumstellar disc lifetime of 1-10 Myr 
which is of the same order to (if not smaller than) the giant 
planet formation timescale. To circumvent this problem. 
Boss (2002) has proposed that giant planets form directly 
from the gravitational fragmentation and collapse of a pro- 
toplanetary disc. However, there are still a number of open 
issues about this scenario such as the formation and sur- 
vival of bound structures since most calculations so far 
have used an isothermal equation of state and/or too low 
resolution. Furthermore, while solving the timescale prob- 
lem, it is not clear at all if the peculiar composition and 
structure of Jupiter and Saturn can be explained within 
this model. More work is definitively required to investi- 
gate this scenario further. 

Since the discovery by Mayor & Queloz (1995) of the 
first extrasolar giant planet at short distance to its star, 
there is mounting evidence that planets might have formed 
at locations which do not necessarily correspond to those 
where they are observed today. Gravitational interactions 
between the growing planet and the gaseous disc lead to 
angular momentum transfer resulting in a net inward mi- 
gration of the planet and possibly gap formation (Lin & 
Papaloizou 1986, Lin et al 1996, Ward 1997, Tanaka et 
al. 2002). For low mass planets, the migration rate is lin- 
ear with mass (type I migration. Ward 1997) . Higher mass 
planets open a gap and the migration rate is set by the vis- 
cosity independently of planetary mass (type H migration, 
Ward 1997). While the general physical understanding of 



the origin of migration is clear, the actual migration rates 
obtained, especially for type I migration, are so short that 
they are inconsistent with the number of extrasolar plan- 
ets actually detected. Type II migration timescales are 
found to lie between 0.1 to 10 Myr, a timescale compa- 
rable to (if not shorter than) the typical lifetime of the 
disc as well as the planet formation timescale in the core- 
accretion scenario. 

Since all relevant timescales (planet formation, disc 
evolution, and migration) are of the same order of magni- 
tude, it appears difficult to obtain a self-consistent picture 
while omitting anyone of these processes. For this pur- 
pose, we have developed a new code, structured around 
three modules. These modules calculate the disc structure 
and evolution (using the simple formalism of a viscosity) , 
including planet migration, the interaction of planetesi- 
mals with the planet atmosphere, and the planet structure 
and evolution. In section 2 we will describe these modules 
in more details and present the numerical tests we have 
performed to validate them. Section 3 will be devoted to 
the effect of migration and disc evolution on formation 
timescales, and to formation models of extrasolar plan- 
ets. Finally, we will discuss these results and conclude in 
section 4. 

2. Equations and assumptions 

Our code to compute the formation of giant planets 
consists of three modules. The first module follows the 
method given by Papaloizou & Terquem (1999, referred 
to PT99 in the following) to calculate the disc structure 
and its time evolution. The second one calculates the in- 
teraction between the infalling planetesimals and the at- 
mosphere of the growing planet; the last one is a plane- 
tary structure and evolution code, written especially for 
this project. We now describe these modules, and present 
some of the tests we have performed. 

2.1. Disc structure and evolution 
2.1.1. Vertical structure 

The first module aims at determining the structure (both 
vertical and radial) of a protoplanetary disc. The method 
is identical to the one used by PT99 and therefore we only 
briefly recall the main points. Cylindrical symmetry is as- 
sumed and therefore the cylindrical coordinates (r, z, 9) 
are a natural choice. The disc is assumed to be thin. For 
each distance r to the star, the vertical structure is calcu- 
lated by solving the equation of hydrostatic equilibrium, 



where z is the vertical coordinate, p is the density and P 
is the pressiu-e. The disc is assumed to be Keplerian with 
fi^ = GMj,/r^, G being the gravitational constant and M* 
the mass of the central star, assumed to be equal to the 
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solar one. This equation is solved together with the energy 
equation: 



dF 9 2 
— = -pM , 
oz 4 



(2) 



which states that the energy production by the viscosity 
V is removed by the radiative flux F, and the diffusion 
equation for the radiative flux: 



F 



3kp dz ' 

where T is the temperature, k is the opacity, and a 
is the Stefan-Boltzmann constant. The viscosity is cal- 
culated using the standard Shakura & Sunyaev (1973) 
a— parametrization v = aC^ /fl where the speed of sound 
Cj is determined from the equation of state. 

The boundary conditions for this part of the calcula- 
tion are the same as in PT99, formally. 



T{z = H) =T{T,^,n,r,M,ua), 



P{z = H) 



F{z = H) = ^Mstn^ 

OTT 

and 

F{z = 0) = 0. 



k(T{z ^ H),P{z = H))' 
3 



(4) 
(5) 

(6) 
(7) 



These conditions depend on three parameters: Tab the 
optical depth between the surface of the disc (z — H) 
and inflnity, Tb the background temperature, and Mst 
the equilibrium accretion rate deflned by Mst = Stti'I] 
where E = pdz is the usual surface density, and 
i> = J^^ vpdz/Yi. The values for Tab and Tb are the same 
as in PT99 (namely 10~^ and 10 K); the steady-state ac- 
cretion rate is a free parameter. As shown in PT99, the 
structure obtained hardly varies with the first two param- 
eters. 

This system of 3 equations with 4 boundary conditions 
has in general no solution, except for a certain value of 
H . This value is found iteratively: equations ^ to O are 
numerically integrated from z = iJ to z = 0, using a 
fifth-order Runge-Kutta method with adaptive step length 
(Press et al. 1992) until F(z = 0) = to a given accuracy. 

Using this procedure, we calculate, for each distance 
to the star r and each value of the equilibrium accretion 
rate Mst, the distributions of pressure, temperature and 
density T(z; r, Mst), P{z; r, Mst), p{z; r, Mst)- 

Using these distributions, we finally calculate the mid- 
plane temperature {Tmid) and pressure {Pmid), as well 
as the effective viscosity v{r,Mst), the disc density scale 
height H{r,Mst) defined by p{z ^ H) ^ e-^/^p{z = 0). 
The surface density 'S{r,Mst) is also given as a function 
of Mst (for each radius). By inverting this former relation, 
we finally obtain relations T„iid{f','^^), Pmidii",^), i^ir,T,) 
and H{r, E) for each value of r (and each value of the 
other parameters a. Tab and Tb). 



2.1.2. Evolution of the surface density 

The time evolution of the disc is governed by the well- 
known diffusion equation (Lynden-Bell & Pringle 1974): 



dE 



3d_ 

r dr 



or 



l§-^rJir)) 



(8) 



where J(r) = -^^(i/Er^/^) is the mass flux (integrated 
over the vertical axis z). This equation is modified to take 
into account the momentum transfer between the planet 
and the disc, as well as the effect of photo-evaporation and 
accretion onto the planet: 



dE 



3d_ 

r dr 



r^l^^yY.r^l'^ 
dr 



Mr) 



-E^(r) + Qpianct(?')-(9) 



The rate of momentum transfer A between the planet and 
the disc is calculated using the formula derived by Lin & 
Papaloizou (1986): 



A(r) = i^^GMst^r 
2r 



Mplanct \ 



AT, 



max r 



IH) 



,(10) 



where a is the sun-planet distance and /a is a numeri- 
cal constant The photo-evaporation term E^, is given by 
(Veras & Armitage 2004): 



= 

cx 



for R < Rg 
ior R> R 



(11) 



9' 



where Rg is usually taken to 5 AU, and the total mass loss 
due to photo-evaporation is a free parameter. Finally, a 
sink term Qpianet is included in Eq. |51 to take into account 
the amount of gas accreted by the planet. This term is 
generally negligible compared to the other ones, except 
during the runaway phases. 

To solve the diffusion Eq. El we need to specify two 
boundary conditions. The first one is given at the outer 
radius of the disc (in our simulations this radius is usu- 
ally taken at 50 AU). At this radius, one can either give 
the surface density E or its temporal derivative. Since the 
characteristic evolution time of the disc is the diffusion 
timescale 



Tj/ cx — cx 

iy a 



(12) 



which^ is proportional to r^/^ for discs of approximately 
constant aspect ratio (which is the case in these models 
(see PT99)) the outer boundary condition has little influ- 
ence. 

The second condition is specifled at the inner radius 
where we have used the following condition: 



(13) 





= 0. 




dr 


inner radius 



^ In this formula, the disc scale height H is the scale height 
of the unperturbed disc, and not the scale height in the middle 
of the gap. 

^ The second part of Eq. 1121 is obtained by expressing the 
equation^as ^ ^ Q,^H and then replacing the sound velocity 
by QH in the definition of v. 
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Since the total mass flux through a cyhndcr of radius r is 
given by: 



<i>(r) = 2TTrJ(r) = Sni'T, + Gnr- 



dr ' 



the boundary condition Eg. 1131 can be expressed as: 



$(r) 



inner radius 



(14) 



(15) 



i. e. the mass flux through the inner radius is equal to the 
equilibrium flux. Therefore, this condition is equivalent to 
say that the inner disc instantaneously adapt itself to the 
conditions given by the outer disc. As discussed in PT99, 
this is consistent with the expression of the characteristic 
timescale as a function of the radius fEa. I12|l . 

2.2. Migration rate 

Dynamical tidal interactions of the growing protoplanet 
with the disc lead to two phenomena: inward migration 
and gap formation (Lin & Papaloizou 1979, Ward 1997, 
Tanaka et al. 2002). For low mass planets, the tidal inter- 
action is linear, and migration is of type I (Ward 1997), 
whereas higher mass planets open a gap, leading to a re- 
duction of the inward migration (referred to type II mi- 
gration). 

Analytical models of type I migration have been com- 
puted by Ward (1997). The resulting migration timescales 
are much shorter than both the disc lifetime and the 
planet growth timescale, making survival of forming plan- 
ets unlikely: the planet is accreted onto the central star. 
Migration could be stopped if there is an inner cavity in 
the disc, but planets at larger distances remain difficult 
to explain. Tanaka et al. (2002) have performed new an- 
alytical calculations of type I migration, in two or three 
dimensional discs and found longer migration timescales 
but still too short to ensure survival. Their migration rate 
is nevertheless confirmed by recent three dimensional nu- 
merical calculations of disc structure and planet migration 
(Bate et al. 2003). 

On the other hand, suggestions of increased type I mi- 
gration timescales can be found in Nelson & Papaloizou 
(2004). As shown by these authors, the torques exerted 
on at least low mass planets (A/pianct < 30 Af®) embedded 
in turbulent MHD discs are strongly fluctuating, result- 
ing in a slowing down of the net inward motion. Contrary 
to laminar discs (as considered by Tanaka et al. 2002 and 
Bate et al. 2003) the migration proceeds as a random walk, 
and the mean value of the migration velocity seems to be 
highly reduced, compared to the laminar case. Moreover, 
as shown by Menou & Goodman (2004), type I migration 
of low-mass planets can be slowed down by nearly one 
order of magnitude in regions of opacity transitions. 

These considerations seem to indicate that the actual 
type I migration timescale may in fact be considerably 
longer than the one originally estimated by Ward (1997) 
or even by Tanaka et al. (2002). For these reasons, and for 
lack of better knowledge, we actually use for type I migra- 
tion the formula derived by Tanaka et al. (2002) reduced 



by an arbitrary numerical factor // chosen between 1/10 
and 1/100. Tests have shown that provided this factor is 
small enough to allow planet survival, its exact value does 
not change the formation timescale but just the extent of 
the migration (see section rrT|l . 

The migration velocity for low mass planets is taken 
to be: 



da 



planet 



dt 



= — 2//apia 



^planet 



(16) 



where ipianot = Mpianct(GM,apia„ct)^^^ is the angular 
momentum of the planet and the total torque F is given 
by: 

F = (1.364 + 0.541as,p) f %^^V Spr|,17^, (17) 

\ JW* Os,P / 

where Cg is the sound velocity and as = ^^7- In this 
expression, the subscript P refers to quantities at the lo- 
cation of the planet. 

For type II migration, two cases have to be considered. 
For low mass planets (when their mass is negligible com- 
pared to the one of the disc) the inward velocity is given by 
the viscosity of the disc. As the mass of the planet grows 
and becomes comparable to the one of the disc, migration 
slows down and eventually stops. In this latter case, the 
variation of the planet orbital angular momentum is equal 
to the angular momentum transport rate in the gaseous 
disc (Lin et al. 1996, Ida & Lin 2004): 



^ [Afpianotapianct^^] = ^'^i'^T 



(18) 



In all cases of type II migration, the migration rate is 
limited by the viscous transport in the disc: 



da 



planet 



dt 



2a 



planet 



Min(l,2Sa2j^_^^^/Mpianet). (19) 



Migration type switches from type I to type II when the 
planet becomes massive enough to open a gap in the disc. 
It happens when the Hills radius of the planet becomes 
greater than the density scale height H of the disc. 

2.3. The planetesimals 

2.3.1. Surface density and physical properties 

The initial amount of heavy elements in the disc is a poorly 
constrained quantity. For this reason, the dust-to-gas ratio 
is varied in our simulations, and takes two values depend- 
ing on the mid-plane temperature of the disc: /q for 
temperatures below 150 K and 1/4/^i/g for higher tem- 
peratures. In principle, the position of the iceline should 
evolve due to the viscous evolution of the disc. However, 
since our treatment of the planetesimals disc is very sim- 
ple, we do not take into account this evolution. 

We assume that due to the scattering effect of the 
planet, the surface density of planetesimals is constant 
within the current feeding zone but decreases with time 
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proportionally to the mass aecreted (and/or ejected from 
the disc) by the planet. The feeding zone is assumed to 
extend to a distance of 4 i?H on each side of the planetary 



orbit, where i?H 



3M. 



1/3 



'^planet 



is the Hills radius 



of the planet. For the inclinations and eccentricities of the 
planetesimals, we use the following prescription (P96): 



/2GM, 



planetesimal 



'^planet y ^planetesimal 



(20) 



where Mpianetesimai and Tpianetesimai are the mass and ra- 
dius of planetesimals, at the location of the planet, and 



e = max(2i, 2 ). 

^planet 



(21) 



Finally, we also take into account the ejection of plan- 
etesimals due to the planet, using the ejection rate given 
by Ida & Lin (2004): 



accretion rate / K;sc.disk \ 
ejection rate V Kiurf .planet / 



(22) 



where T4sc,disk = ^2GM0/apianot is the escape veloc- 
ity form the central star, at the location of the planet, 
Kiurf, planet = ^GMpianct/^^c is the planet 's characteristic 
surface speed, and Rc is the planet's capture radius (see 
next section). 

It is worth noting that our model for the evolution 
of the disc of planetesimals remains a very simple one 
in which a number of effects are neglected. For example, 
wc! omit the effect of gas drag on the planetesimals which 
given their assumed size of 100 km is reasonable. Using the 
gas density obtained in our model, we can calculate the 
typical timescale for radial drift. For regions above ~ 4 
AU, we find values of order of ~ 10^ to 10^ the orbital 
time. In regions closer to the sun and at the begining of 
our simulation, this timescale may be much lower owing to 
the higher gas density. However, these inner regions evolve 
rapidly to lower densities, and gas drag becomes quickly 
negligible. 

We also neglect, apart for ejection and accretion onto 
the planet, the perturbations created by the growing 
planet in the planetesimal disc such as heating, gap for- 
mation and shepherding (Tanaka & Ida 1997, Tanaka & 
Ida 1999, Thommes et al. 2003). For the latter effect, we 
have calculated the critical migration rate, below which 
shepherding occurs (Tanaka & Ida 1999) and found that 
in all cases the migration rate exceeds this critical value. 
The planet acts as a "predator" and not as a "shepherd" 
(see Tanaka & Ida 1999). 

Apart from their radius, the planetesimals are charac- 
terized by a number of bulk properties like their density 
Pb, their tensile strength or, heat of ablation (vaporiza- 
tion or melting) Qati and some more material dependent 
parameters. This allows the simulation of the fate of differ- 
ent planetesimal types in the atmosphere of the growing 
protoplanet, as described in the next part. 



2.3.2. Interaction with the growing atmosphere 

Given a core and the structure of the surrounding atmo- 
spheric envelope (pressure P, temperature T, density p 
etc. as a function of the distance from the planetary cen- 
ter), the second module computes the trajectory and de- 
struction of planetesimals inside this region by integrating 
a system of ordinary differential equations using a fifth- 
order Runge-Kutta method with automatic step size con- 
trol (Press et al. 1992). This approach is similar to the one 
described by Podolak et al. (1988, hereafter referred to as 
PPR88). We will present this module and its results in 
more details in a oncoming separate paper (see Mordasini 
et al. 2005), and restrict ourselves here to a description of 
the physical aspects we take into account. There are four 
main mechanisms controlling the fate of a planetesimal: 

Gravity: The gravitational attraction is calculated in 
the two-body approximation (planet - planetesimal) as- 
suming a spherical mass distribution. This is justified be- 
cause inside Rn, where the calculation takes place, the 
effect of the planet generally largely predominates other 
the third body effect. 

Aerodynamic drag force: Apart from gravity, the aero- 
dynamic drag force Fd is the second force that defines 
the trajectory of the planetesimal. Using the standard for- 
mula for the drag force at high Reynolds numbers (see e.g. 
Landau & Lifshitz 1959) we have: 



Fd = ^CdPV^S, 



(23) 



where v is the velocity, S is the instantaneous cross sec- 
tion of the planetesimal and Cd is the drag coefficient 
(for a sphere), computed as a function of the local Mach 
and Reynolds numbers using the equations of Henderson 
(1976). These equations give smooth and continuous val- 
ues for Cd for widely varying aerodynamic environments 
(free molecular flow - continuum flow; hypersonic flow - 
incompressible flow). 

Thermal mass loss: The basic mechanism of thermal 
ablation is very simple: the drag force leads to a dissipation 
of kinetic energy, some fraction of which is used to heat the 
gas, while the remainder goes into heating up the body. 
When the heat flux is sufficiently high, the surface of the 
body becomes so hot that melting or vaporization starts. 
From these considerations, the simple classical ablation 
equation known from the study of terrestrial meteors is 
(Opik 1958, Bronsthen 1983): 



dM 



IchPv'S-^. 

^ Qabl 



(24) 



Qabi is the amount of energy per unit mass needed to 
bring body material from its initial temperature to the 
point where melting or vaporization occurs plus the spe- 
cific heat needed for this phase change. Ch, the heat trans- 
fer coefficient, is an unknown function which depends on 
the velocity of the particle, the flow regime, the shape of 
the body etc. and describes the fraction of the incoming 
kinetic energy flux of the gas that is available for ablation. 
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One of the problem with the classical equation - apart 
from the neglect of re-radiation from the hot surface and 
heat conduction into the interior - is that the rate of ther- 
mal mass loss is heavily dependent upon the value of Ch 
which can vary by several orders of magnitude depending 
upon the size of the infalling body (from 10~^ to w 0.5, 
Svetsov et al. 1995). 

For terrestrial, cm-sized meteorites a mean heat trans- 
fer coefficient of about 0.1 can reasonably explain the ob- 
servations (Bronsthen 1983). Whether this value can be 
extrapolated to impactors of tens or even hundreds of kilo- 
meters in size remained for a long time unclear (PPR88). 
Fortunately, the large number of studies of the impact 
of comet Shoemaker-Levy 9 (SL9) onto Jupiter have now 
shown that modelling such an impact with the above equa- 
tion setting Ch » 0.1 greatly overestimates the mass loss 
and predicts too small a penetration depth in compari- 
son to detailed hydrodynamic simulations (Sekanina 1993, 
Field & Ferrara 1994, Ahrens et al. 1994). One can con- 
clude from the SL9 event that for large, km-sized bod- 
ies, thermal ablation is of minor importance compared to 
mechanical destruction by aerodynamical forces (Svetsov 
1995, see below) and that Ch is small for such large ob- 
jects, of the order of 10"'*— lO^"* (sec e.g. Field and Ferrara 
1995, Zahnle & Mac Low 1994). Anyway, it is clear that 
thermal ablation must be considered in all cases even when 
fragmentation occurs, as it controls the final dissolution of 
the impactor (Korycansky et al. 2000), especially if one is 
interested not only in the energy deposition profiles (which 
was of major interest in the SL9 case) but also in the mass 
deposition (as it is the case here, since we aim to calculate 
the fraction of mass dissolved inside the atmosphere, and 
the actual mass of the solid core). 

In order to obtain a realistic model of the energy in- 
put, we follow Zahnle (1992) by replacing the conductive 
heat influx term per surface unit a \pv^ in the supersonic 
contimmm flow by C h .rad{Tps)<jT^g whenever the latter is 
smaller, a way to take into account that in this regime, 
the energy input is proportional to the radiative flux from 
the shock front (Opik 1958). In this equation, Tps is the 
post-shock temperature which is directly computed as a 
function of the local atmospheric properties and the im- 
pactor velocity by solving numerically the normal shock 
wave jump conditions for a real gas (Landau & Lifshitz 
1959), using the equation of state of Chabrier et al. (1992). 
Our results for Tps arc very similar to the ones obtained 
by Chevalier et Sarazin (1994). CH,radiTps) is a coefficient 
which denominates the fraction (< 1) of radiation reach- 
ing the planctcsimal surface due the screening effect of 
ablated material. We use here the results of Bibermann et 
al. (1980). 

We also include, where appropriate, radiation from the 
undisturbed atmosphere (PPR88) and a convective heat 
input mechanism which is proportional to the temperature 
difference between the surrounding gas and the impactor 
surface temperature. We use the expressions presented by 
Sibulkin (1952), corrected -if necessary- for a turbulent 



boundary layer. A similar expression can also be found in 
Svetsov (1995). 

The actual surface temperature of the planctcsimal 
and the amount of ablated mass are obtained by equat- 
ing the energy influx to the amount of re-radiation and 
energy used for ablation. For the reasons presented in 
PPR88 we neglect the heat conduction towards the inte- 
rior of the planetesimal. For vaporization, dM/ dt is linked 
to the wall temperature via the Knudsen-Langmuir for- 
mula (Bronsthen 1983, PPR88). 

The thermal mass loss rates found for large, non- 
fragmented bolides are in our simulation considerably 
lower than in a model using the simple ablation equation 
and Ch = 0.1. 

Mechanical destruction: In deeper layers of the at- 
mosphere, the aerodynamic pressure acting on the front 
side of the impactor can overcome the internal strength 
of the body (either tensile strength or self-gravity, see 
PPR88). The bolide will then approximately act as a fluid 
and undergo deformation, fragmentation and mechanical 
ablation (Svetsov et al. 1995). Numerical simulations of 
the SL9 event have proven the prevalent importance of 
these effects to understand the atmospheric destruction 
of a large impactor (see e.g. Zahnle & Mac Low 1994). 
These effects occur in two different regimes. In the first 
one (static regime, when the aerodynamic pressure load- 
ing builds up on a time-scale larger than the time required 
for a sound wave to travel through the body, see Svetsov 
et al. 1995, Korykansky et al. 2000), we use a model which 
combines the effects of lateral spreading (using the so 
called "pancake" model of Zahnle 1992) and the growth 
of Ray leigh- Taylor instabilities (Sharp 1984, Roulsten & 
Ahrens 1997, Korycansky et al. 2000) at the front side 
of the fluidized impactor to obtain a multi-staged frag- 
mentation model. In the second one (dynamical regime, 
when the pressure loading is dynamical) , where no lateral 
spreading is observed (Svetsov et al. 1995, Korycansky et 
al. 2000) but where instead mass is removed at the front 
of the bolide, we use a simple model for mechanical abla- 
tion which is also based on the growth of RT instabilities. 
For large impactors which are found in this regime, mass 
loss by mechanical ablation overcomes thermal ablation 
usually by a large factor. 

After the calculation, the three dimensional energy and 
mass deposition profiles are converted back into one di- 
mension to make them usable by our third module de- 
scribed in the next section. The planetesimal impact code 
is also used to calculate the capture radius, a calculation 
performed in a same way as described in P96. 

2.4. Protoplanet structure and evolution 

2.4.1. Equations of the internal structure 

The third module calculates the planet internal structure 
including a growing core (inciluciing the energy deposited 
by accreted planetesimals) and the gas accretion due to 
both the contraction of the envelope and the increase of 
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the outer radius of the planet. The standard equations of 
planet evolution are solved: 



dm 
dL 
dm 
dP 

dm 
dT 

IP 



3 

47rp' 
dU P dp 



dt 



dt 



Vad or Vrad, 



(25) 
(26) 
(27) 
(28) 



using opacities from Bell & Lin (1994) and Alexander & 
Fergusson (1994), and the equation of state (EOS) from 
Chabrier et al. (1992). In these equations, r, L, P, T are re- 
spectively the radius, the luminosity, the pressure and the 
temperature inside the atmosphere. These four quantities 
depend on the gas mass m included in the sphere of ra- 
dius r. p is the mass density, U the specific internal energy, 
and Mcore the mass of the solid core, eacc is the amount 
of energy released by the accretion of planetesimals (see 
section l2.r>.2|l : this term dominates largely the energy bud- 
get in Eq. 1261 until the maximum accretion rate has been 
attained (see Section 12. 4. The temperature gradient is 
given by the adiabatic (Vad) or by the radiative gradient 
(Vrad), depending on the stability of the zone against con- 
vection which we check using the Schwarzschild criterion. 
For convective zones, we assume that the temperature gra- 
dient is given by the adiabatic one; in other words, we do 
not use the mixing length theory (MLT) in these models. 
Test have shown that including MLT does not change our 
models. 

Two approximations can be used regarding the fate 
of planetesimals' matter, after their destruction inside the 
envelope. In the "sinking" approximation (see P96), this 
mass is assumed to slowly sink toward the core, leading to 
an extra term in the core luminosity. In the "no sinking" 
approximation, the matter is assumed to remain inside the 
envelope. Note that, as in P96, the matter is added to the 
core mass in Eq. |^in both cases. Under this aspect, the 
sinking approximation appears to be more self-consistent. 

Two inner boundary conditions are necessary to 
solve Eq. [53 to the core radius is set to i?coro = 

1/3 



and 



Tsmi — 7mid(5j(aplanct5 ^))j 



(30) 



where Pmid(S(apianot, t)) and rmid(S(apianet, i)) are cal cu- 
lated using the structure of the disc (see section 12. 1.1(1 . 

This condition is valid as far as the disc can supply 
enough mass to keep the outer radius equal to the Hill (or 
the accretion) radius, i.e. if the gas accretion rate is below 
the maximum accretion rate (see section 12. 4. 2|l . 

2.4.2. Solid accretion rate 

We use the expression of Greenzweig & Lissauer (1992) to 
calculate the solid accretion rate: 



dM 



solids 



dt 



(31) 



where Sp is the surface density of solids at r = flpianot and 
Rc is the capture radius of the planet. The capture radius 
is calculated using our second module (see section l2.3.2f) . 



2.4.3. Gas accretion rate 

Normally, the gas accretion rate is determined from the 
condition: 



R{m = Af, 



planet } 



min(i?H, R a 



(32) 



where i?accr is the accretion radius of the planet (see P96). 
At each timestep, the mass of the envelope (and then the 
total mass) is determined by iterations to satisfy this con- 
dition to a given accuracy. From comparing models at t 
and t + dt, we obtain the gas accretion rate. 

The condition R{m = Mpianct) = nrin(_RH, -Racer) can 
be fulfilled only if the disc can supply enough gas to the 
planet. Once a gap opens in the disc, the maximum gas 
accretion rate is set to the rate given by (Veras & Armitage 
2004): 



gas, max 



= 1.668 (Mpianct/Mj)'/'e ^^+0.04, (33) 



disc 



(iT^Sr) ^^'^ *li*' ^'^^^ luminosity Lcorc is given by where Mj is the Jupiter mass and M, 



the sum of the remaining kinetic plus the correspond- 
ing potential energy of planetesimals, after having passed 
through the atmosphere. The density of the core is fixed 
to 3.2 g/cm^ as in P96. 

The mass of the envelope is given by requiring that 
the outer radius of the planet is equal to the minimum of 
the Hill radius and the accretion radius i?accr = "^^'^Sf""* 
where Cg is the sound velocity inside the disc at the loca- 
tion of the planet (see P96). 

Finally, two outer boundary conditions are necessary. 
They are given by requiring that the disc and the planet 
join smoothly at the outer radius: 



Pmid (S(apianct , t)) 



(29) 



is the disc ac- 
cretion rate away from the planet. When the maximum 
accretion rate is reached, the growth of the planet mass is 
set by the disc and no longer by the internal structure of 
the planet, which is no longer computed. 



2.5. Initial conditions and physical assumptions 

We start our calculations with a core of O.6M0, at a dis- 
tance astait from the star. The initial surface density of 
the disc is usually taken as a power law, the total mass 
being given (for fixed boundaries) by the normalisation at 
5.2 AU. Its life time is given by the value of a and the rate 
of photo-evaporation. At last, the solid surface density is 
fully characterized by the dust to gas ratio Jd/g- 
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Table 1. Physical parameters used in our models 



Parameter 



Value 



Core density 
Initial core mass 
Helium fraction 
Internal disc radius 
External disc radius 
Planetesimal radius 



3.2 g/cm^ 
0.6M© 

0.24 
0.25 AU 
50 AU 
100 km 



In addition, two numerical parameters have to be spec- 
ified, // (reduction of type I migration) and /a, the nu- 
merical factor in the expression of the momentum trans- 
fer between the planet and the disc. As shown in Alibert 
et al. (2004), modifications in the outer planet bound- 
ary conditions due to the gap formation have a small ef- 
fect. Therefore, the gap essentially limits the gas accretion 
onto the planet, and we will restrict ourselves to the case 
/a = 0. 

Some other parameters of the model, which are not 
changed in the calculations presented here are summarized 
in Table n 

2.6. Numerical tests 

Various tests have been performed to validate our entire 
model. Concerning the disc part of our code, we have com- 
pared the resulting functions i>(r, S) for different values of 
a to the analytical fits given in PT99. Figure ^ shows a 
comparison for a = 10"^. We recall that PT99 mention a 
difference between the numerical results and their fits of 
the order of 50% at most. This is also comparable to the 
maximum difference between our results and the fits. 

The interior structure module has been used to repro- 
duce the results of BP86. Figure |21 shows the evolution of 
a forming Jupiter, for a planet below the cross-over mass 
(until the mass of accreted planetesimals and the mass of 
accreted gas are equal). In this calculation, the accretion 
rate of planetesimals is constant, and there is no inter- 
action between the envelope and the planetesimals. The 
boundary conditions are those of case 1, 6 and 7 of BP86. 
Figure |21 shows the central temperature and density for 
these models, which is in good agreement with their re- 
sults (see Figure 1 in BP86). We conclude that our internal 
structure module works properly. 

Our planetesimal accretion module has been tested ex- 
tensively by comparing results of simulations of impacts 
into terrestrial and Venusian atmosphere (ReVelle 1978; 
Hills & GodaHnH Zahnle 1992), as weU as into Jupiter's 
atmosphere (Zahnle & Mac Low 1994). Figure |3| shows 
one of these tests: energy deposition profiles are plotted 
for compact ice comets of different initial radii hitting the 
Jovian atmosphere. The initial conditions are chosen to 
match the SL9 event. As expected, larger bolides pene- 
trate deeper and have a larger peak energy deposition. 
Our results are virtually identical to those of the analytical 
model of Zahnle & Mac Low (1994), which is not surpris- 



00 
o 




Fig. 1. Mst(r, S) given by our calculations (solid lines), 
and the analytical fits by PT99 (dotted lines), for a = 
10~^. The values of the radii (in AU) at which these rates 
are computed are indicated on the plot. 
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Fig. 2. Central density and temperature for models analog 
to case 1, 6 and 7 of BP86. 



ing as the pancake model is used in both cases. The edge 
at an altitude of about 90 km comes from the radiative 
limit for thermal ablation and appears at the same point 
as in Zahnle & Mac Low (1994) and Crawford (1996). The 
thick line schematically represents an energy release profile 
obtained from a high-resolution hydrodynamic simulation 
by Mac Low & Zahnle (1994), predicting a very similar 
peak energy deposition altitude and value. A model using 
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Fig. 3. Energy deposition profiles of SL9-type impactors 
with different initial radii in the Jovian atmosphere. For 
comparison, a schematic representation of the energy de- 
position profile of a 0.5 km impactor in the high-resolution 
hydrodynamic simulations of Mac Low & Zahnle (1994) 
is given (thick line). 



the classical ablation equation and a high heat transfer 
coefficient {Ch ^ 0.1), but no mechanical effects can not 
reproduce this profile. 

Finally, the entire code has been tested using the same 
initial conditions as P96 (case Jl) turning disc evolution 
and migration off. In this case, we obtain a cross-over time 
of ~ 8 Myr, in close agreement with their result (compare 
Fig. 0] with Fig. 1 of P96). We conclude that our code 
properly reproduces giant planet formation in absence of 
migration and disc evolution. 

3. Results 

3.1. Formation timescales 

As shown in BP86 and P96, the formation timescale is es- 
sentially given by the time needed to reach the cross-over 
mass. To quantify the effects of migration and disc evo- 
lution on planet formation, we compare different models, 
all reaching the cross-over phase at the same location in 
the disc, namely 5.5 AU. 

For comparison purpose with P96, we consider an ini- 
tial disc density profile given by S cx r^^, like in P96, again 
for comparison purpose. The constant is chosen to yield 
E = 525 g/crn^ at 5.2 AU; this corresponds to ~ 3 times 
the surface density in the minimum mass solar nebula at 
the position of Jupiter. This surface density profile yields 
isolation masses that do not depend on the distance to the 
sun. We do not take into account photo-evaporation, and 
we start with an embryo of 0.6 initially at 7, 8 or 15 
AU, depending upon the choice of the parameter // (0.01, 
0.03 and 0.1 respectively). The viscosity parameter a is 
set to 2.10^^ and the dust-to-gas ratio is equal to 1/70 
for disc mid-plane temperature below 150 K, and 1/280 




2x10" 4x10^ 6x106 
time (yr) 



xlO'' 



Fig. 4. Mass of accreted planetesimals (Afhoavy, solid line), 
mass of accreted gas (dotted line), and total mass of the 
planet (dashed line) as a function of time, for conditions 
similar to case Jl of P96. 



in the opposite case. Finally, we use for these simulations 
the no-sinking approximation (see P96). 

Figure [51 shows the mass of accreted planetesimals and 
function of time for the models considered here. 
Note that the mass of accreted planetesimals does not cor- 
respond to the core mass since some fraction of them are 
being destroyed while traversing the envelope and never 
reach the core. For the in situ case (without disc evolu- 
tion), the time to reach the cross-over mass is around 30 
Myr, much longer than the typical discs lifetimes (Haisch 
et al. 2001). For the simulations with migration and disc 
evolution, the corresponding timescales are respectively 
~ 1, ~ 1 and ~ 2.2 Myr for an embryo starting at 7, 8 or 
15 AU respectively. 

Allowing for migration and disc evolution, we obtain a 
time to reach the cross-over mass of about 1-3 Myr, i.e. 
more than ten times faster than in our identical model 
in which migration and disc evolution have been switched 
off. As explained in Alibert et al. (2004), the main reason 
for this speed-up is found to be due to migration which 
prevents the severe depletion of the feeding zone as ob- 
served in the in situ formation model. Migration therefore 
suppresses the need to accrete a gas envelope in order to 
extend the feeding zone (phase 2 of P96) and the planet 
reaches cross-over mass much faster. 

As stated in section 12.3.21 the effect of ablation is 
found to be negligible, disruption of planetesimals occur- 
ring very deep in the envelope. The two cases, sinking and 
no-sinking, considered by P96 differ therefore much less in 
our calculations than in P96. For the simulation started 
at 8 AU, the time to reach the cross over mass is found to 
be of order ~ 1.2 Myr using the sinking approximation, 
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Fig. 5. Mass of accreted planetesimals (lines starting at 
O.6M0), and mass of accreted gas as a function of time, 
until the cross-over mass is reached. Solid lines: starting 
at 8 AU (with // = 0.03), dotted lines: starting at 15 AU 
(with // = 0.1) and dashed lines: in situ model. The model 
starting at 7 AU is very close to the one starting at 8 AU 
and is not shown. 

compared to ~ 1 Myr in the no-sinking case. For that rea- 
son, and since it is more self-consistent, we will use in the 
following the sinking approximation, as in P96's standard 
calculations. 

3.2. Formation of extrasolar planets 

Depending on the initial parameters used, very different 
planets can form. Figure [3 shows the results of a set of 
1000 simulations'^, in a mass vs semi-major axis diagram. 
The initial parameters for these simulations are Ostart — 
3, 4, 5, 6, 7.5 10, 15 and 20 AU, fo/c = 1/25, 1/50, 
1/70, 1/100 and 1/200, and a photoevaporation rate of 
lO-^Mo/yr, 2x lO-^Mg/yr, 4x lO-^Mg/yr, lO-^Mg/yr 
and 2 x lO~®M0/yr. The initial disc is given by S cx r~'^/^, 
the normalization at 5.2 AU being between 100, 300, 500, 
700 and 900 g/cm^. This corresponds to disc masses of 
0.01 to ~ O.IMq between 0.25 and 50 AU. The viscosity 
parameter is fixed to 2 x 10^"^, and we use // = 0.03 for the 
whole set of simulations. The edges of the disc are at 0.25 
and 50 AU, so planets whose inner radius of the feeding 
zone is below 0.25 AU are considered to have migrated 
inside the star, and are represented by dots on the vertical 
axis. 

One must keep in mind that the diagram shown here 
does not take into account any statistical weighting of 
the different initial parameters, so one should be careful 

^ Each simulation takes ~ 3 to 4 hours on a modern PC. 



Fig. 6. Distance to the central star for the same models 
as in figure [S] The kinks around 0.8 and 2 Myr signal 
the change from type I to type II migrations. Solid line: 
starting at 8 AU (with // = 0.03), dotted line: starting at 
15 AU (with // = 0.1), and dashed line: in situ model. 

when comparing the obtained distribution with the ob- 
served one. We present it simply to illustrate the possible 
diversity of planets that can form in our model. In a forth- 
coming paper we will discuss the probability of occurrence 
of these planets. 

As an illustration taken from this set of simulations, 
we provide some details on the formation process of 
one particular planet. This object evolved from the fol- 
lowing: Ostart = 15 AU, a disc photo-evaporation rate 
4 X 10~^Mo/yr, an initial normalisation of the disc at 
5.2 AU equal to 500 g/cm^, and fo/c ^ V^O. With a fi- 
nal mass of ^ 3.5Mj, and a final distance ~ 2.5 AU from 
its star, this planet is indicated by the circle in Figure 
Figure IHI shows the evolution of the mass of accreted gas, 
the mass of accreted planetesimals, as well as the mass 
of the disc, and Figure |51 gives the distance to the central 
star as a function of time. The cross-over mass is reached 
after ^1.6 Myr, and shortly after, due to gap formation, 
the accretion rate of gas is limited to its maximum value, 
which decreases with decreasing disc mass. The formation 
process ends after 5.5 Myr when the disc has essentially 
disappeared. Note however, that 90% of final mass of the 
planet is reached already after ~ 4 Myr. At the end of 
the simulation, the planet as accreted ~ 30Mq of plan- 
etesimals. As stated in Section 13.11 this mass does not 
correspond to the mass of the core, since some of the ac- 
creted planetesimals may have been destroyed during the 
crossing of the envelope. Our model provides an estimate 
of the core mass by tracking the location where the plan- 
etesimals are destroyed (second module). For this simu- 
lation, ~ 6.4M© of accreted planetesimals reach the core 
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Fig. 7. Final mass and semi- major axis (dots), for a set of 
simulations, using different initial parameters (see text for 
details). The crosses are the parameters of the observed 
exoplanets as in November 2003. The big circle represents 
the final parameters of the model detailed in section 
The dashed line gives the inner limit of our simulations 
(when the inner radius of the planet feeding zone is below 
0.25 AU). 

without being destroyed in the envelope. Ignoring further 
processes, this provides the mass of the core which is sim- 
ilar to the one inferred from Jupiter's internal structure 
models (Guillot et al. 2004). The total content in heavy 
elements will finally depend on the metallicity of the ac- 
creted gas. For solar composition (prior to the formation of 
planet esimals) , the accreted gas add an additional ^ 7Mq 
of heavy elements. 

The migration of the planet can be divided into three 
phases. Before ~ 1 Myr the planet undergoes type I mi- 
gration at which time a gap opens and migration switches 
to type II. Shortly after ^ 2 Myr, the mass of the planet 
becomes non negligible compared with the disc mass and 
migration slows down and eventually stops at the time the 
disc disappears. 

4. Summary and discussion 

We have presented a new code devoted to the calcula- 
tion of giant planet formation models, taking into account 
the protoplanetary disc structure, as well as the migra- 
tion of the forming planet. These calculations show that 
the formation of giant planets, at least in the first phase 
until runaway gas accretion, can be heavily sped-up if one 
takes into account the effect of migration. This is mainly 
due to the suppression of phase 2 as described in P96. 
Using an initial disc model similar to the one of P96, we 
obtain a time to reach the cross-over mass of the order of 



Fig. 8. Mass of the different components for model of sec- 
tion l3.2l Solid line: mass of accreted planetesimals, dashed 
line: mass of H/He, and heavy solid line: mass of the disc. 




2xl0« 4x10*^ 6xl0« 

time (yr) 



Fig. 9. Distance to the sun for model of section IX^ The 
kink around 1 Myr corresponds to the change from type I 
to type II migration. 

1 Myr, significantly shorter than the typical disc lifetimes. 
Moreover, this speed-up due to migration has been found 
to be robust against reasonable changes in various param- 
eters (Alibert et al. 2004). Therefore, migration not only 
accounts for the orbital distribution of extra-solar planets, 
but also considerably shortens the formation timcscale in 
the core-accretion scenario. The formation of giant plan- 
ets in this scenario is therefore in excellent agreement with 
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inferred disc lifetimes without having to consider discs sig- 
nificantly more massive than the minimum mass solar neb- 
ula. 

We note that the speed-up due to migration does not 
preclude other effects to further reduce this timescale. For 
example, reducing the dust opacity would decrease the 
critical mass (Ikoma et al. 2001), thus leading to another 
reduction in the formation timescale. Furthermore, such 
a reduced opacity could account for the existence of giant 
planets with small central core. 

Using different initial conditions can lead to the for- 
mation of a wide variety of planets. However, the compar- 
ison of our results with observations of extrasolar plan- 
ets need to take into account the probability distribution 
of various initial conditions. Work is under progress to 
obtain, using a Monte-Carlo approach, synthetic distribu- 
tions which can then be compared in a statistical way with 
the observed distributions. 

These calculations are however subject to some un- 
certainties, among them the simplified treatment of the 
planetesimals disc, or the calculation of the ejection rate, 
that directly determines the final heavy elements content. 
Work is under progress to improve these aspects. 

Finally, it seems to be very difficult to form a planet, 
and to prevent it from spiraling into the sun if the amount 
of type I migration as computed today is not reduced by 
a factor of at least 10. 
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